########################################################
## PROGRAM NAME: 004_dmv_map.R                        ##
## AUTHOR: MATT MLECZKO                               ##
## DATE CREATED: 04/26/2021                           ##
## INPUTS:                                            ##
##    002_af_wide_final.Rda                           ##
##                                                    ##
##                                                    ##
## OUTPUTS:                                           ##
##                                                    ##
## PURPOSE: Create Figures 1 and 2                    ##
##                                                    ##
## LIST OF UPDATES:                                   ##
########################################################

log <- file("004_dmv_map.txt") 
sink(log, append=TRUE)
sink(log, append=TRUE, type="message")

## load libraries ##
library(tigris)
library(tidyverse)
library(viridis)
library(sf)

## define paths ##
data_path <- "PATH TO DATA HERE"
progs <- "PATH TO PROGRAMS HERE"

## set working directory ##
setwd(data_path)

## load analytic file ##
load("002_af_wide_final.Rda")

##
## create map of integrated vs. segregated tracts in DC metro area ##
##

## get shapefiles ##
dc2010.geotrts <- tracts(
  state = "DC",
  #county = NULL,
  year = 2010)

va2010.geotrts <- tracts(
  state = "VA",
  county = c("059","153","107",
             "013","510","179",
             "177","061","047",
             "683","187","630",
             "610","043","113",
             "157","685","600"),
  year = 2010)

md2010.geotrts <- tracts(
  state = "MD",
  county = c("031","033","021",
             "017","009"),
  year = 2010)

wv2010.geotrts <- tracts(
  state = "WV",
  county = c("037"),
  year = 2010)

##
## bring on data from analytic file ##
##

incoming <- af.wide.final %>%
  select(GEOID,
         trt_entropy_ers_2020,
         trt_entropy_sess_2020,
         trt_entropy_joints_2020,
         ld_er_2020,
         md_er_2020,
         hd_er_2020,
         ld_ses_2020,
         md_ses_2020,
         hd_ses_2020,
         ld_jt_2020,
         md_jt_2020,
         hd_jt_2020,
         int_er_relf_2020,
         int_ses_rel_a_2020,
         int_jt_relf_a_2020)

summary(incoming$trt_entropy_joints_2020)

##
## combine all the data ##
##

dmv2010.geotrts <- rbind(dc2010.geotrts,
                         va2010.geotrts,
                         md2010.geotrts,
                         wv2010.geotrts) %>%
  rename(GEOID = GEOID10) %>%
  left_join(incoming, "GEOID") %>%
  mutate(FIPS = paste0(STATEFP10,
                       COUNTYFP10)) 

summary(dmv2010.geotrts$trt_entropy_joints_2020)

##
## main plots ## 
##

##
## figure 1 ##
##

ggplot() + 
  geom_sf(data = dmv2010.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## figure 2 ##
##

ggplot() + 
geom_sf(data = dmv2010.geotrts, 
        aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## now just PG county ##
##

pg.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "24" & 
         COUNTYFP10 == "033")

ggplot() + 
  geom_sf(data = pg.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = pg.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## now just Montgomery county ##
##

mont.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "24" & 
           COUNTYFP10 == "031")

ggplot() + 
  geom_sf(data = mont.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = mont.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

#
## now just Fairfax county ##
##

ff.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "51" & 
         COUNTYFP10 == "059")

ggplot() + 
  geom_sf(data = ff.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = ff.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## now just Loudoun county ##
##

lou.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "51" & 
           COUNTYFP10 == "107")

ggplot() + 
  geom_sf(data = lou.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = lou.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## now just Arlington county ##
##

arl.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "51" & 
           COUNTYFP10 == "013")

ggplot() + 
  geom_sf(data = arl.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = arl.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## now just Alexandria county ##
##

alx.geotrts <- dmv2010.geotrts %>%
  filter(STATEFP10 == "51" & 
           COUNTYFP10 == "510")

ggplot() + 
  geom_sf(data = alx.geotrts, 
          aes(fill = trt_entropy_joints_2020)) + 
  scale_fill_gradientn(colors = plasma(5),
                       name="Joint diversity",
                       na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

ggplot() + 
  geom_sf(data = alx.geotrts, 
          aes(fill = as.factor(int_jt_relf_a_2020))) +
  scale_fill_manual(name = "Joint integration",
                    breaks = c("0","1"), 
                    labels = c("Segregated","Integrated"),
                    values=c("purple4", "gold"),
                    na.value = "grey100") +
  theme_minimal() + 
  coord_sf(datum = NA) 

##
## alternative plots ## 
##

counties <- counties(cb = TRUE,
                     year = 2010)

dmv.counties <- counties %>%
  mutate(FIPS = paste0(STATE,
                       COUNTY)) %>%
  filter(FIPS %in% dmv2010.geotrts$FIPS)

ggplot(dmv.counties) +
  geom_sf(aes(fill = CENSUSAREA)) +
  geom_sf_label(aes(label = NAME))

##
## more plots ## 
##

plot(dmv2010.geotrts[c("trt_entropy_ers_2020",
                       "trt_entropy_sess_2020",
                       "trt_entropy_joints_2020")])


plot(dmv2010.geotrts["int_er_relf_2020"])
plot(dmv2010.geotrts["int_ses_rel_a_2020"])
plot(dmv2010.geotrts["int_jt_relf_a_2020"])

plot(dmv2010.geotrts[c("int_er_relf_2020",
                       "int_ses_rel_a_2020",
                       "int_jt_relf_a_2020")])


### END OF PROGRAM ###

sink()
sink(type = "message")

